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Abstract 



In this article we study the numerical approximation of incompressible miscible displace- 
ment problems with a linearised Crank-Nicolson time discretisation, combined with a mixed 



finite element and discontinuous Galerkin method. At the heart of the analysis is the proof 



of convergence under low regularity requirements. Numerical experiments demonstrate that 
the proposed method exhibits second-order convergence for smooth and robustness for rough 
problems. 



\q ■ 1 Introduction and Initial Boundary Value Problem 

in 

Mathematical models which describe the miscible displacement of fluids are of particular economical 
relevance in the recovery of oil in underground reservoirs by fluids which mix with oil. They also 
play a significant role in CO2 stratification. 

This publication extends the analysis of [1] , which studies the discretisation of miscible displacement 
under low regularity. Unlike to [1] which is based on a first-order implicit Euler time-step (leading 
to a nonlinear system of equations in each time step), here we examine the discretisation in time 
by a linearised second-order Crank-Nicolson scheme. Crucially, the new, more efficient method 
inherits stability under low regularity. Like in [1], the concentration equation is approximated with a 
discontinuous Galerkin method, while Darcy's law and the incompressibility condition is formulated 
as a mixed method. High-order time-stepping for miscible displacement under low regularity has 
recently also been addressed in [4] , however, with a continuous Galerkin discretisation in space and 
discontinuous Galerkin in time. We refer for an outline of the general literature to [1, 2, 3, 4]. 

Definition 1 (Weak Formulation). A triple (u,p,c) in 

L°°(0,T;ifjv(div;ft)) X L°°(0, T; L%(Q)) x (L 2 (0, T; H 1 ^)) n H 1 ^, T; H 2 ({1)*)) 
is called weak solution of the incompressible miscible flow problem if 
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(Wl) for t € (0, T), v £ H N (div; Q) and q € L§(«) 

(//(c) K _ V v) - (p,dwv) = (p{c)g,v) 

(g,divu) = (q 1 - q P ,q). 

(W2) forallw£@{0,T;H 2 (n)) 
T 

- (<j) c, d t w) + (D(u) Vc, Vw) + (it • Vc, w) + (/c, w) - (cq 1 , w) dt = 0. 

(W3) c(0,-) =c in H 2 (n)*. 

For the data qualification we refer to condition (A1)-(A8) in [1] and for the physical interpretation 
of the system to [1, 2, 3]. We point out that D growths proportionally with u: 

d {l + \u\)\i\ 2 < t T B(u,x)t< d°(l + |n|)|e| 2 , u,£ € R d , x£fl. 

Thus D is in general unbounded on Lipschitz domains f2 and in the presence of discontinuous 
coefficients, which are permitted in this paper. 

2 The Finite Element Method 

We compactly recall the definition of the finite element spaces from [1]. Let = to < t\ < . . . < 
tM = T be a partition of the time interval [0, T]. Let kj := tj — tj-\ and dta- 1 := k~ x [a? — a^ 1 ). 
We consider meshes 7 of 17 with elements K and set hx '■= diam(K). We denote by S S (T) the 
space of elementwise polynomial functions of total or partial degree s. For Wh € S S (T) the function 
VhWh is defined through {VhWh)\K = ^{wh\k)- The sets of interior and boundary faces are 
£q(T) and 8,gfi(7). We set £(T) = £q(T) U £an(T) and assign to each E G £(T) its diameter 
He- We denote jump and the average operators by [•] and {•}. The concentration c is discretised 
at time j on the mesh 7 3 c or simply by 7 3 . The approximation space for the variable c at time 
step j is denoted by $ 3 C . Often we abbreviate £ J := £(Tc), £^ := £n(Tc), £qq := £-dd(7i). We 
denote the Raviart-Thomas space of order £ by RT^(Tu). The approximation spaces of u and p 
are §1 := RT^(Ot) n H N (div;Q) and S£ := S*(o£) n Lg(fi). We frequently use the global mesh 
size and time step /i- 7 := max^ gT j uT j /i/^, /i := maxo<j<M , A; := maxo<j<M fe J as well as to 

S n = n^li^M'Sp = Iljli §P)§c = lljlo^c- In addition we impose conditions (M1)-(M5) of [1] 
which are on shape-regularity, boundedness of the polynomial degree, control Hi^Hl 4 ^5 ll^/illi? 1 an d 
the structure of hanging nodes. 

To deal with discontinuous coefficients and the time derivative, we substitute D by 

D h : L 2 (n) d -> § s (7 c ,R dxd ), v ^ n T oD(ti,-) 

where the Ilg- are projections such that 1 1 II 2D 1 1 j<- < Given quantities a J , a J_1 and a- 7 " 2 at 

times tj, tj-i, tj-2, we denote a? = \a? + \a?~ l and a = fa- 7 " 1 — \a?~ 2 . 
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The diffusion term of the concentration equation is discretised by the symmetric interior penalty 
discontinuous Galerkin method: Given Ch,Wh G §c, Uh G Si, we set 

B d {c h ,w h ;u h ) := (B^u/JV^, V h w h ) - ([c h ], {D J h (u h ) V h w h }) ^ 

- (Hd\{ D £("/0 V^c^}) £ ^ + (<r 2 [c ft ], [Wh])^ 

where a is chosen sufficiently large to ensure coercivity of B^, cf. [1]. The convection, injection and 
production terms are represented by 

B cq (c h ,w h ;u h ) := 1 / 2 {{u h V h c h ,w h ) - {u h c hl V h w h ) + ((q 1 + q p )c h ,w h ) (1) 

where (uh ■ n) + : = max{ii^ • n, 0} and (uh • n)- : = min{u^ • n, 0}. We set B = B d + B cg . 
Algorithm (A dG ). Choose 4 t e §i for j = 0, 1. Given find (*4,_p{) G §>i x § J P such that 



(M4) K lu h,v h ) ~ {Ph,divv h )= (p(4) 9,v h ), 
(q h ,dwu h ) = ((q 1 - q P y,qh)- 



(2) 



For 2 < j < M find c* h G Sc swc/i i/mi, /or a// u>h £ Si, 

(4>d t c j h ,w h ) + B{ck> ,w h ;u{) = (^q I3 ,w h ) (3) 
anc? so/ue to obtain (u J h ,p J h ) G Si x Sp. 

The algorithm only requires the solution of a linear system in each time step. The iterate c\ can 
be computed with an implicit Euler method and fine time steps. The use of extrapolated values 
such as u J h is classical, e.g. see [5, p. 218]. 



3 Unconditional Well-posedness, Boundedness and Convergence 

Given c^ _1 and cj^ 2 , there exists a solution G Si of (3) because the bilinear form B is positive 
definite. For t G [tj-i, tj], let Ch(t, •) := c? h + tj ^- c? h x . Then d t Ch{t, •) = d t c' h (-). We interpret 

elements of § u , § p and S c as time-dependent functions with stepwise constant values. Let 

\ c h\l h ■= (Phiuh^hCh^hCh) + {o- 2 [ch], i c h\) E J n + {\uh ■ n&j\ [c h ], [ch])^- 

Theorem 1. Let p° = ||p||oo- There exists a constant C > such that 

\K\\ + ||divfi» || + \\j? h \\ < (|| A|| + M 1 - q P \\) (4) 

holds for all j = 2, 3 . . . , M. Equally we have 

ii^vji 2 + f my dt < ii^/2 c i|| 2 + /■* w^Y^fdt ( 5 ) 

for all j = 2,3..., M. 
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Proof. The stability of v? 1 , v? 2 , p 3 l ,p* 2 follows from a classical inf-sup argument. This implies 
stability of W and p> . We choose Wh = cK 1 in (3) to verify that 

dtU^clf + \c£\\ + W + q P ) l,2 cZ\\ 2 < 2(4>d t clcf?) + 2B(^,cT;^) = 2(?q I \cR i ). 
The Cauchy-Schwarz inequality, multiplication by ki and summation over i give 

II^Mf + E^t* < ll^ 1/2 4ll 2 + £M(^ 1/2 ^ll 2 

for all j = 2, 3, . . . , M. □ 

For simplicity the next theorem is stated assuming meshes are not adapted in time. For the 
extension to changing meshes consult [1]. However, observe that that the discretisation with the 
implicit Euler method gives additional stability in kiW^^dtc^W 2 , which allows to change meshes 
more rapidly. 

Theorem 2. The time derivative dfCh belongs to L 2 (t\,T; H 2 (VL)*) and 

l|5tC/ l || i 2( iliT . H 2(Q)*) = \\d t C h \\ L 2( tuT . H 2( n )*) < 1, 

independently of the mesh size and time step. 
Proof. Let Wh € $ J C . We recall from [1] 

Ba^w^ul) < (1 + \\ui\\ 1/2 ) \<? h y (||V h ti;ft|| L 4( n ) + \\w h \\ L i {n) + IHttfjJIIgi ), 
B cq (c i h ,w h ;u J h ) < (1 + \\u J h \\ h ) \<? h y (\\VhW h \\ + ||i«h||£4(n) + \\a[w h ]\\^), 

lkK]|| 2 4 <(i + KII)^ 1/2 lkll^ ( n)- 

With L 2 -orthogonality and 



/ ((pd t c j h ,w)dt= (c 3 q I3 ,w h ) - B(c h J ,w h ;u J h ) dt 

Jti Jt x 

< J o (l + Kll)(l + ll<lli(div;Q))l4Lill^ll^(n)di 



~ \\ w \\L 2 (0,T;H 2 (n)) 



one completes the proof. □ 

Theorem 3. Let (ui,pi, Ci)i^ be a sequence of numerical solutions with (hi, hi) — ► as i — > 00. 
Then there exists c € L 2 (0,T;H 1 (n))nH 1 (0,T;H 2 (n)*) such that, after passing to a subsequence, 
d^c inL 2 (Q T ), d t Ci d t c inL 2 (0,T;H 2 (Q)*) andVa Vc in L 2 (0,T; H^(n)). lfc\,c\ -» c 
m H 2 (Q)* then c satisfies (W3). 
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The proof is, up to the treatment of the initial conditions, exactly as in [1]. It is based on the 
Aubin-Lions theorem and the embedding 



§, s (7i) [BV(fi) n L\Q), L\Q)} 1/2 L 2 (ft), 
where [•, denotes the complex method of interpolation. 

Theorem 4. Lei (ui,pi, Cj)j g ]N fre numerical solutions with (hi,ki) — > and q — > c m L 2 (S7r) as 
? — > oo. There exists u G L°°(0, T; i?Ar(div; 0)) and p G L°°(0, T; Lq(O)) suca taat, a/ter passing to 
a subsequence, Ui —> u in .ffjv(div; fi) and pi — > p in Lq(Q) as (hi,ki) — > 0. Furthermore, (u,p,c) 
satisfies (Wl). 

Proof. Use Strang's lemma, for details sec [1]. □ 
We interpret Ui as piecewise constant function in time, attaining in (tj-i,tj] the value |n(t J_1 ) — 



Theorem 5. Let (ui,pi,Ci)i e ^ be a sequence of numerical solutions with (hi,ki) — > as i — > oo 
and letu G L°°(0, T; Hn(<Hv ; fi)) and c G L 2 (0, T; tf 1 ^)) n ^(O, T; H 2 (Q)*) be a limit offacfii 
in the sense of Theorems 3 and 4- Then (u,c) satisfies (W2). 

Proof. Let v G f^(0, T; and t>j(t) G §c an approximation to ?;(i) in (tj-i, tj]. Using the 

strong convergence of (VhVi)i in L°°(Q,T) d and the weak convergence of the lifted gradient of q in 
L 2 {tt T ) d , we find 



I (Vc,D{u)Vv)dt = lim / (V h Q, D /l (« i )V h «*) - ([q], {B h (ui)V h Vi}) £n dt. 



As in [1] it follows that Bd(ci,Vi;ui) coincides in the limit with (Vc, B(u)Vd). One can also 
conclude by adapting [1] that 

/ (u • Vc, u) + (a 7 c, v) dt = lim / B cq (a, vf, Uj) dt. 
7ti i-KWti 



One arrives at 



/ -(<f)C,d t v) + (B(u)Vc, V?;) + (u- Vc,?;) + (o 7 c,?;) - (cg 7 ,?;)dt 
lim / {(t>d t c 3 h ,w h )+B{c^,w h \u 3 h )-{^q I \w h )dt = {). 



'*1 

Hence (W2) is satisfied for u G 9(0, T; if 00 (ft)). The extension to £F(0, T; H 2 (Q)) follows from 
boundedness and density of smooth functions. □ 
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K = k 2 Id 


K = k[ Id 


IK = k, Id 

* 






Figure 1: Example 1: Left: computational domain; right: absolute value \uh\ of the Darcy velocity 
at t = 1.0 before any interaction between the concentration front and the corner singularity. 





Figure 2: Snapshots of Ch at t = 1.5 and 2.0, computed with the Crank-Nicolson scheme. 
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4 Numerical Experiments 



The numerical experiments are carried out in two space dimensions with the lowest-order method 
on a mesh which consists of shape-regular triangles without hanging nodes and which is not changed 
over time. The diffusion-dispersion tensor takes the form 

D(«, x) = <f>{x) (d m U + |u| d t E(u) + |u| d t (Id - E{u))) . (6) 




Numerical Example 1 (Singular Velocities). To examine the effect of a singular velocity field 
caused by a discontinuous permeability distribution and a re-entrant corner we employ the L- 
shaped domain and K with k\ = 0.1 and &2 = 10 -6 as depicted in Figure 1. The injection 
and production wells are located at (1,1) and (0,0), respectively. The porous medium is almost 
impenetrable in the upper left quarter, forcing a high fluid velocity at the reentrant corner where 
the nearly impenetrable barrier is thinnest. This leads to a singularity |u| ~ r~ a , where r is the 
distance to the reentrant corner and a ~ 1, cf. [1]. Figure 2 shows the concentration when the 
front passes the corner and at a later time. The solution c\ t contains steep fronts but shows only 
the localised oscillations that are characteristic for dG methods. 

Numerical Example 2 (Convergence rates). Convergence rates are determined by comparing 
the numerical solution to a reference solution c re f that is computed with high accuracy on a one 
dimensional grid. More precisely, we set (f) = 1, c = 1, K = l and g = and choose f2 to be the 
ball .6(0, 1) C K 2 . Using polar coordinates (r, 93), we choose q 1 = 4(1 — r) 6 and q p = |r 6 . Then 
the Darcy velocity only changes in the radial direction and is determined by an ODE, which has 
the nonnegative exact solution u(r) = £ (3r 6 — 24 r 5 + 70 r 4 — 112 r 3 + 105 r 2 — 56 r + 14). Conse- 
quently, the concentration equation reduces to a linear parabolic equation in one space dimension. 
Figure 3 shows snapshots of the solution c re f with d m = 1.0 x 10~ 5 , dt = 4.0 x 10~ 4 and Figure 
4 shows that 1? error of implicit Euler method is of order 0(h 2 + k) whereas the Crank-Nicolson 
reaches the order 0(h 2 + k 2 ). 
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